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Abstract 

One of the more promising recent approaches to turbulence modelhng is the Varia- 
tional Multiscale Large Eddy Simulation (VMS LES) method proposed by Hughes 
et al. [Comp. Visual. Sci., vol. 3, pp. 47-59, 2000]. This method avoids several con- 
ceptual issues of traditional filter-based LES by employing a priori scale partitioning 
in the discretization of the Navier-Stokes equations. 

Most applications of VMS LES reported to date have been based on hierarchical 
bases, in particular global spectral methods, in which scale partitioning is straight- 
forward. In the present work we describe the implementation of the methodology 
in a three-dimensional high-order spectral element method with a nodal basis. We 
report results from coarse grid simulations of turbulent channel flow at different 
Reynolds numbers to assess the performance of the model. 
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1 Introduction 



Large-eddy simulations (LES) provides a physically more appealing framework 
for turbulent flow prediction than the more traditional Reynolds-averaged 
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models (RANS). In the latter the full impact of the ensemble averaged effect 
of turbulent advcction on the mean flow field has to be modeled. The essence 
of the LES approach on the other hand is to directly solve (with a complete 
time and space resolution) the three-dimensional and time-dependent motion 
of the largest turbulent scales. These scales are in general associated with the 
most energetic motion of the turbulence field and it is (ideally) only the least 
energetic motion that need to be modeled. The concept as such is therefore 
well suited to confront the scale complexity and transient behavior inherent to 
turbulent flows and offers a more complete representation than RANS models 
per se. 

In traditional LES, large- and small-scale motion are separated by applying a 
spatial filtering operation to the Navier-Stokes equations. This results in a set 
of equations for the large-scale motion. The residual motion, i.e. turbulent mo- 
tions on scales that are smaller than the filter width, appear in these equations 
as a residual stress term that must be modeled. There arc several conceptual 
issues in filter-based LES that have to be addressed. For instance, filtering 
and spatial differentiation do not in general commute on bounded domains or 
for non-uniform grids, and it is not obvious how to prescribe correct bound- 
ary conditions for the filtered velocity at solid walls. Another unwarranted 
character of filter-based LES models is that the residual stress model has a 
tendency to affect the entire range of the spectrum and not only represent the 
filtered effect of the unresolved scales near the spectral cut-off. These issues 
have been the subject of a considerable amount of research, and the lesson 
learned, in general, is that LES works well in cases where the rate-controlling 
processes occur at the largest (resolved) scales of motion, or equivalently in 
flows where the unresolved scales, and consequently the model, only plays a 
secondary dynamical role. 

In this paper we consider a different approach to LES, the variational multi- 
scale (VMS) LES method originally proposed by Hughes et al. [1]. The VMS 
LES method employs an a priori scale partitioning in the discretization of 
the Navier-Stokes equations, instead of flltering to separate the large- and 
small-scale motion. The scale partitioning appears to overcome some of the 
disadvantages of fllter-based LES. First, since there is no flltering, all issues 
concerning commutation errors and boundary conditions at solid walls are ad- 
dressed. Second, since the scale partitioning is performed during discretization, 
we develop different equations representing different ranges of the spectrum. 
Different modelling assumptions can then be applied to each range of the spec- 
trum, improving our ability to apply the model terms where they are needed, 
and only there. 

We implement the VMS LES formulation in a high order spectral element 
method for the solution of the Navier-Stokes equations. Spectral element meth- 
ods offer an attractive combination of the accuracy of spectral methods and 
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the flexibility of finite element methods. This provides us with an attractive 
framework for model development in which the numerical errors can be con- 
trolled, such that the true performance of the model can be assessed. The first 
implementations of the variational multiscale LES method [2,3,4] used global 
spectral methods. These methods naturally employ an orthogonal modal basis, 
such that the scale partitioning becomes straightforward. Recently, the method 
has also been implemented in the context of other numerical schemes, such 
as finite element methods [5,6] and finite volume methods [7]. Our spectral 
element code uses an element-wise discretization with nodal basis functions 
that contain information on all the scales. One of the challenges of the present 
work is therefore to devise a way to separate the large and the small scales, 
and to implement the VMS terms. We show that this can be achieved by an 
element-by-element transformation into the Legendre modal basis functions. 

In the following sections we will discuss the variational multiscale method as a 
turbulence modelling tool, and describe the implementation of the method in a 
spectral element solver for the incompressible Navier-Stokes equations. Finally 
we will present computed results, from both a high-resolution DNS and coarse 
grid VMS LES for the turbulent flow in a plane channel at different Reynolds 
numbers. The computed results show that, even with simple modelling applied 
to the small-scale equations, the performance of the methodology is promising. 



2 The Vciriational multiscale method 



In this section we will discuss the variational multiscale method as a tool for 
turbulence modelling. The variational multiscale LES method was introduced 
by Hughes et al. [1] and later elaborated by Collis [8]. We will outline the 
method following Collis, to shed light on the modelling assumptions employed 
in the derivation of the model. 

The Navier-Stokes equations describing the dynamics of a viscous, incompress- 
ible fluid are 

V • w = 0, (la) 

— + w • Vu = -Vp + uV^u + /, (lb) 
at 

where the independent variables are the velocity, u— {u,v,w), and the pres- 
sure, p. The kinematic viscosity is denoted by i/, and / is a body force term. 
The non-dimensional parameter that characterizes the flow is the Reynolds 
number Re = |u|L/i/. 

For ease of presentation we assume homogeneous Dirichlet boundary condi- 
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tions for the velocity, i.e. 

u{x) = X er. 

We can then construct the weak, or variational, formulation by choosing test 
and trial functions in the same function space V. Note however that in general 
the test and trial spaces will differ at the boundary. 

U = {u,p) eV 
W = {w,q) eV 

We take the inner product of W with Eq. (1) (written in the compact form 
M{U) — F) to obtain the weak Navier-Stokes operator: 

{wMm = C{W, U) - n{w, u) = {W, F), (2) 

comprising the linear Stokes operator 

du 

C{W, U) = {w, - (V • w,p) + {V'w, 2i/V^u) + (r, V • u), (3) 
and nonlinear advection represented by the Reynolds projection 

TZ{w,u) = B{w,u,u), (4) 
where B is the tri-linear term. 

B{w,u,v) = {Vw,uv). (5) 

To take into account the multiscale representation, we write the solution space 
V as a disjoint sum 

v = vevev, 

in which V and V comprise the large and small scales, respectively, whereas V 
contains the unresolved scales that cannot be represented by the numerical dis- 
cretization. The scale partitioning is sketched in Fig. 1. Now, by decomposing 
the test and trial functions in these spaces 

u = u + u + u, 
w + w + w, 

we can develop exact variational equations governing different scales. Further- 
more, by assuming that the scale partitioning is orthogonal, we obtain the 
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Fig. 1. Schematic of the turbulent energy spectrum with scale partitioning 
following equations governing the large, the small, and the unresolved scales 

U) - n{w, u) - (W, F) - n{w, u) - C{w, u, u) 



= TZ{w, u) + C{w, u, u) + C{w, u, u), 

C{W, U) - 7^(^^, u) - {W , F) - 7^(^^, u) - C(n, n, u) 
= TZ{w,u) + C{u,u,u) + C{w,u,u), 

C{W, U) - n{w, u) - C{w, u, u) - C{w, u, u) 

= n{w,u) +n{w,u) +C{w,u,u) + iW,F), 



(6a) 
(6b) 
(6c) 



where C{w, u, u') = B{w, u, u') + B{w, u', u) is the cross stress term. We 
have written these equations in a form such that all terms that depend on the 
unresolved scales are collected in the right-hand sides. It is thus evident that 
there is an effect of the unresolved scales on the computable, resolved scales, 
and it goes without saying that this effect must be modeled. In the original 
paper by Hughes et al. [1], the modelling assumptions were not stated, but 
the issue was clarified by Colhs [8], who showed that essentially the following 
assumptions result in a method that is identical to the method proposed by 
Hughes (which by then had produced excellent results [2,3]) 

• The separation between large and unresolved scales is sufficiently large so 
that there is negligible direct dynamic influence from the unresolved scales 
on the large scales. 

• The dynamic impact of the unresolved scales on the small scales are on 
average dissipative in nature. 



The simple scalar Smagorinsky-type model is in an averaged sense fully consis- 
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tent with the last assumption. In order to approximate the temporal behaviour 
at the cut-off, a more refined modelling approach would be needed. This is 
however outside the scope of the present study. 

With these assumptions, the LES model is only applied to the small scale 

equation, adding additional dissipation where it is mostly needed. Different 
implementations of this method by the Hughes group [2,3], by Ramakrishnan 
and CoUis [4], and by Jeanmart and Winckelmans [9] have produced very good 
results even for wall-bounded channel flows. 

We remark here that both assumptions are, or at least should be, open to 
scrutiny. Firstly, although it is plausible that the unresolved scales do not 
influence the large scales, it is not necessarily obvious. In fact, a recent anal- 
ysis by Oberai et al. [10] showed that the energy transfer from the large and 
small scales, respectively, to the unresolved scales depends critically on the 
discretization method and the function spaces that are used to perform the 
scale partitioning. Furthermore, Reynolds number effects or other aspects of 
the flow physics may mandate that a more sophisticated model for the large 
scales must be taken into account. Secondly, the assumption of a one-way 
cascade from the small to the unresolved scales require that flow is properly 
resolved, such that the cut-off is far out in the inertial range. This is unfor- 
tunately not always the case in LES computations. Such considerations are, 
however, outside the scope of the present study. Our objective is to present an 
implementation of the VMS LES formulation in the spectral element method. 
For this purpose the assumptions employed to date [1,8] are acceptable. At 
present, we merely note in passing that the VMS method presents an excellent 
framework for improved modelling to address these issues. 

Bearing the above in mind, we can formulate the variational modeled equa- 
tions. The effect of the unresolved scales on the large scales, given by the 
right-hand side of (6a), is neglected according to the first assumption, while 
the effect of the unresolved scales on the small scales, given by the right- 
hand side of (6b), is modeled by a Smagorinsky term. The equation for the 
unresolved scales is naturally omitted. The resulting set of equations is 

£(TF, U) - n{w, u) - n{w, u) - C{w, u, u) - {W, F) = 0, (7a) 

C{W,U) -n{iv,u) -n(w,u) -C(w,u,i^^ , , 

~ (7b) 
- {yV, F) = -(V'mT, 2utV'u). 

The terms that couple the different scales are evident in (7); the small-scale 
equation has been supplemented with a dissipative term that accounts for the 
interactions between the small and the unresolved scales, whereas large-scale 
Reynolds and cross stress projection account for the large-scale influence on 
the small scales. The large-scale equation contains a projection of the small- 
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scale Reynolds stress onto the large-scale to account for interaction between 
the small and the large scales (i.e. back-scatter). 

We are however chiefly concerned with the complete resolved solution U ~ 
U + U, not with the large and small scales per se, and adding the large- and 
small-scale equations we obtain 

(w, J^(U)^ + {V'w, 2i/r V"w) = (W, F) . (8) 

We note that in this equation, all the interaction terms between the large 
and the small scales arc accounted for in the advection operator TZ, which is 
part of the first term on the left-hand side in (8). The projected cross and 
Reynolds stress terms that appear in the large- and small-scale equations (7) 
are therefore mainly important for analysis and turbulence modelling, but need 
not necessarily impact on the implementation of the method. The variational 
formulation is hence primarily an analysis tool and a vehicle for developing the 
VMS methodology. The essential feature of the method is that the turbulence 
modelling should be confined to the small scales. As long as a suitable scale 
partitioning can be performed on the solution space, the methodology can in 
principle be applied to any discretization, as indicated by Hughes et al. [2]. 



3 Implementation in the spectral element method 

In this section we describe the implementation of a VMS LES model in a high- 
order spectral element method for the solution of the incompressible Navier- 
Stokes equations. 

We will start with a brief discussion of Legendre polynomials and the spectral 
element basis functions. These concepts are important, both for the descrip- 
tion of the basic numerical method as well as for the implementation of the 
variational multiscale framework that follows. More details about the topics 
covered in Sections 3.1 and 3.2 can be found in [11]. 



3.1 Legendre spectral elements 

The Legendre polynomials are orthogonal with respect to the unweighted inner 
product in the function space L^(— 1, 1). The Legendre polynomials are given 
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by the recurrence relation 



Lo{x) = 1, 

Ll{x) = X, /gX 

2k + 1 k 

Lk+i{x) = ———xLk{x) - -——Lk-i{x), k>l, 
k + 1 k + 1 

where Ln{x) is the Legendre polynomial of degree A^". 

The Gauss-Lobatto-Legendre (GLL) points {■Cil^o ^ = 1] defined 
as the extrema of the A^th order Legendre polynomial Ln{x), in addition to 
the endpoints of A: 

Co = -1, Cat - 1, L'^iQ = 0, J = 1, . . . , AT - 1, eo < 6 < • ■ • < Civ. (10) 

Furthermore, the Gauss- Legendre (GL) points {r)j}^~-^ on A, that are used to 
represent the pressure in the spectral element method, are defined implicitly 
by LN-i{Tij) = 0, i.e. as the zeros of the Legendre polynomials of order (A^ — 
1) [12]. Note that the GL points do not include the endpoints of A. 

The spectral element nodal Gauss-Lobatto-Legendre basis is defined by choos- 
ing trial and test functions to be the corresponding Lagrangian interpolants 
at the Gauss-Lobatto-Legendre (GLL) grid points, constructed as A^th order 
polynomials such that each function has the property 

hj{Ci)^Sij, i^O,...,N. (11) 

A function w{x) defined on A can then be represented by the interpolating 
polynomial: 

N 

Wh{3') = ^Wihi{x), X e A, (12) 

i=0 

where Wi = w(^j) are the function values at the GLL points. Higher-dimen- 
sional trial and test functions are constructed as tensor products of these one- 
dimensional functions. Each velocity component is represented this way on 
each element, and the global representation is the sum of the representations 
on all elements. 

A Gauss-Legendre nodal basis for the pressure is constructed in an analogous 
way, only taking into account that we use lower-order polynomials in the basis 
for the pressure to avoid spurious pressure modes in the solution [13] . 



3.2 Spectral element Navier-Stokes solver 



To solve the Navier-Stokes equations (1) we employ an implicit-explicit time 
splitting in which we integrate the advective term explicitly, while we treat 
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the diffusive term, the pressure term, and the divergence equation imphcitly. 
After discretization in time we can write (1) in the form 



{al - i/V^)m"+^ = Vp + g{f, w", u^-\ ...), (13a) 
V • = 0, (13b) 

in which the explicit treatment of the advection term is included in the source 
term g. In the actual implementation we use the BDF2 formula for the tran- 
sient term, 

dt 2At ^ 

which gives a = 3/2At in (13), while we compute the advective contributions 
according to the operator- integration- factor (OIF) method [14] . 

The spatial discretization is based on a spectral element method [13,15]; the 
computational domain is sub-divided into non-overlapping hexahedral cells or 
elements. Within each element, a weak representation of (13) is discretized by 
a Galerkin method in which we choose the test and trial functions from bases 
of polynomial spaces 

u'l e Pn{x) ® PN{y) ® PNiz), (14a) 
p'' e Pn-2{x) ® PN-2{y) ® Pn-2{z). (14b) 

The velocity component variables are defined in the Gauss-Legendre-Lobatto 
basis described above, and they are C°-continuous across element boundaries. 
The pressure variable is represented in the Gauss-Legendre basis, and is dis- 
continuous across element boundaries. As we noted above, the unknowns, or 
dependent variables, in the discrete formulation are the function values of the 
velocities in the GLL points, and of the pressure in the GL points. 

The GLL grid corresponding to the Legendre polynomial of degree N has 
{N + 1) points. Gauss-Lobatto-Legendre quadrature at the {N + 1) GLL 
points is exact for polynomials of degree up to {2N — 1). Hence, the com- 
putation of the inner products corresponding to the diffusive terms in (1) are 
calculated exactly, whereas the evaluation of the non-linear advective terms 
incurs quadrature (aliasing) errors. These errors can be detrimental to the sta- 
bility of the method and must be controlled. The most fundamental approach 
to de-aliasing is to perform over-integration [16,17] - that is, to over-sample 
by a factor 3/2 and calculate the quadrature at this refined grid for the inner 
products containing non-linear terms. The overhead involved depends on the 
amount of the total computational time that is originally spent on the advec- 
tion part, but for the channel flow calculations presented here, over-integration 
typically leads to an increase of around 20% of computational time. 

An alternative, and computationally more efficient approach, is to use polyno- 
mial filtering of the solutions as proposed by Fischer and Mullen [18], where 
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a simple filter operator with negligible computational cost is applied to the 
solution at every time-step. The effect in the spectral space on each element 
is to transfer a certain fraction (the filter strength) of the energy on the high- 
est order basis polynomial in each element over to the third-highest order 
polynomial [19]. By this operation, the pile- up of energy on the highest order 
polynomial is reduced, while the values at the element boundaries are un- 
changed. Filter strengths as small as 1-5% can have positive effects on the 
solution. 

For the solution of the discrete system of equations we now introduce the 
discrete Helmholtz operator, 

where A and B are the global three-dimensional stiffness- and mass matrices; 
the discrete divergence operator, D\ and the discrete gradient operator, G. Ap- 
propriate boundary conditions should be included in these discrete operators. 
This gives the discrete equations 

- = Bf, (15a) 

-L>w"+^ = 0, (15b) 

where the change of sign in the pressure gradient term is caused by an inte- 
gration by parts in the construction of the weak form of the problem. This 
discrete system is solved efficiently by a second order accurate pressure cor- 
rection method. If we let Q denote an approximate inverse to the Helmholtz 
operator, given by a scaled inverse of the diagonal mass matrix, the pressure 
correction method can be written 

Hu* = Bf + Gp"", (16a) 
DQG{p'^+^ - p") = -Du* (16b) 
= w* + gG'(p"+^ - p"), (16c) 

where u* is an auxiliary velocity field that does not satisfy the continuity 
equation, i.e. Du* ^ 0. 

The discrete Helmholtz operator is symmetric and diagonally dominant, since 
the mass matrix of the Legendre discretization is diagonal, and can be ef- 
ficiently solved by the conjugate gradient method with a diagonal (Jacobi) 
preconditioner. Whereas the pressure operator DQG is easily computed; it is 
ill-conditioned. The pressure system is solved by the preconditioned conjugate 
gradient method, with a multilevel overlapping Schwarz preconditioner based 
on linear finite elements [20]. 
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3.3 Incorporation of VMS LES in the SEM 



The implementations of the variational multiscale LES method reported in [2,3,4] 
used global spectral methods. These methods naturally employ an orthogo- 
nal modal basis, such that the scale partitioning becomes straightforward. Our 
spectral element code uses on an elcmcnt-wisc discretization based on the Leg- 
endre polynomials. The Legendre polynomials offer an orthogonal hierarchical 
basis. Like the majority of the spectral element community, we do however 
use a nodal basis constructed from the Lagrangian interpolant functions. In 
this case all the basis functions contains information on all the scales and the 
scale partitioning is no longer straightforward. 



3.3.1 Nodal and modal bases 

We have demonstrated above that in the nodal Gauss- Lobatto- Legendre basis 
a function w{x) defined on — 1 < x < 1 can be represented by a combination 
of the interpolating polynomials, as given by (12). The coefficients in the sum 
are the function values at the grid points. 

An alternative, modal, representation is to use an expansion directly in the 
Legendre polynomials 



'2i + l 



N 

0=0 V ^ 



where the unknowns now are the spectral coefficients Cj. The factor \/^^ 
is used to normalize the basis. The scaled Legendre polynomials represents a 
natural orthonormal basis, in which it is straightforward to perform the scale 
partitioning. In this setting, it is natural to associate the low order polynomials 
with the large scales and the higher order polynomials with the smaller scales. 

The nodal and modal bases are related through the linear transformation 

Kc = w, (18) 

where the entries of the matrix K are given by 



and c and w are the vectors of spectral coefficients and GLL nodal function 
values, respectively. 

Let N = N + N, such that N is the dimension of the polynomial basis for the 
large scales and N is the dimension of the small-scale space. The large-scale 
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Fig. 2. Large- and small-scale partitions in the 2-dimensional polynomial wavenum- 
bcr space. The chosen partition operators arc shown to the right. 

part of a nodal function w can then be written as 

w^KTK-^w, (19) 

where T — diag(/^, 0;^) is the operator that annihilates the small-scale compo- 
nents in the modal basis. For notational convenience, we define the large-scale 
extraction operator 

L = KTK-\ 

while the corresponding small-scale extraction operator is 

S^I-L. 



When tensor products of these operators are formed in higher dimensions, the 
resulting operators extract the components with large-scale, or small-scale, 
respectively, components in all dimensions. The sum of these two operators 
does not add up to the identity, so we choose to define the three-dimensional 
small-scale extraction operator to be 

S = I -{L,®Ly®L^). (20) 

This is illustrated in two dimensions in figure 2. The resulting small-scale 
extraction operator returns functions with small-scale structure in at least 
one dimension. 



3.3.2 Properties of the large- small partition 

The large-scale extraction operator corresponds to a sharp cut-off in the Leg- 
endre modal space. To illustrate the effect in Fourier space of this large-small 
partitioning, we represent the function 

f{x) = f: cos{kx) (21) 

fc=0 
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Fig. 3. Fourier representation of a sharp cut-off in Legcndrc modal space. 

on a spectral element grid on [0, 27r] with 6 elements and 7 grid points in 
each element. Higher wave-numbers can not be represented accurately on this 
particular grid. We extract the large- and small-scale partitions using N — Aoi 
the 7 modes (57%) on each element as the large-scale space. The two resulting 
functions are sampled on a 54-point regular grid, and their Fourier spectra are 
plotted in figure 3. 

The main point illustrated by figure 3 is that although the scale partitioning 
in the Legendre space is done as a sharp cut-off, the Fourier spectra of the 
two partitions are much smoother. The reason for this is that each the original 
cosine terms is represented by a combination of local Legcndrc modes on each 
element. We also note that the gradual growth in the small-scale spectrum 
starts around the cut-off percentage, so the impact of the small-scale extraction 
is weaker than for a straightforward Fourier representation. 

In a more general case with variable element size and/or polynomial order, it 
may be possible to vary the cut-off point in the local Legendre space to keep 
the corresponding "average wavelength" approximately constant throughout 
the whole domain. 



3.3.3 Implementation of the model term 

We now turn our attention to the implementation of the variational multiscale 
model term (V^to, 2i/tV*'S) from (8). Note that the turbulent eddy viscosity 
ut is not a material property of the fiuid, but a property of the fiow field and 
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as such varies through the flow domain. 

It is instructive to first consider the one-dimensional case with a constant 
eddy viscosity. Furthermore, for ease of exposition, we only consider a single 
element. In this case, the weak form of the model term above is 

dw du 



Using the small-scale extraction operator defined above, we have 

N N 

^{^) = YSmqU'^hm{x), (23) 

m=0 q=0 

and for a given test function on Lagrange form {w'^{x) = hi{x)) 

N 

w\x) = Y.Sp,h,{x). (24) 

p=0 

Inserting these representations and using Gauss-Lobatto quadrature, we ob- 
tain 

N N N N 

r=0 p=0 m=0 q=0 
N N N N 

p=0 m=0 q=0 r=0 
N N N ^25) 

— ^ ^ ^pi ^ ^ ^ ] SmqU'^-^pm 
p=0 m=0 q=0 
N 

= J2{S^AS)igu'^ = S^ASu 

q=0 

= (/ - LfAu, 

where the final line is in the form we generalize to higher dimensions. It is easily 
seen from the second-to-last line that (25) represents a symmetric operator 
acting on u. 

The corresponding term in three dimensions is 

(Vw, Vui) = ((5" ® 5^ ® A^) - {L''^ ® ly'^ ® L^^)(5" ® 5^ ® A^)) Ui 

+ ((5^ (g) 5^) - {L'^ ® L^^ ® L'-'^){B' ^A^^ B^)) 

+ [{A' ® (8) S^) - {L'^ (8) L^^ (g) L^^)(A^ (8) (8) S^)) 

(26) 

for each component Ui. 



Ui 

U.: 
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Taking into account that the eddy viscosity, utIx, y), is not constant but rather 
a function that varies in space, will distort the tensor product structure of (26). 
Following the procedure for discretization of terms with variable coefficients 
described in [11], we can write 



Ui 



(Vw, 2i/t(x, y, z)Vui) ^2(l'®P® D^^) V {F ®P® D^) ■ 

- 2 (l^^ ® L^^ (8) (D^L^)^) V {F ®P® D^) Ui 
+ 2 [F® D^^ ® r) V {F ® (8) F) Ui 

- 2 (l"^ ® {DyL^f ® L^^) V {F ® (g) F) Ui 

+ 2 (d"^ ® f ® r) \/ {D' ® F ® r ) 5i 

- 2 ({D'Ff ® L^^ (8) L^^) (L>^ (g) P (8) P) Ui. 

(27) 

In this equation, D denotes the GLL derivation matrix in each direction. Fur- 
thermore, the values of the eddy viscosity are lumped with the GLL integration 
weights in the diagonal matrix V with the entries u^^^p^p^p^, in which rst are 
the grid point indices and V is ordered to be consistent with the ordering of 
the element grid points. 

We are now finally ready to consider the model term in the form given in (8). 
Since the product of a symmetric and an anti-symmetric tensor is zero, we 
find that we only need to to compute the inner product 

In tensor product form, the VMS small-scale dissipation term for the compo- 
nent of the momentum equation becomes 
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(29) 



and we obtain similar expressions for the other two components. The couplings 
between the velocity components, introduced by the second term of (28), are 
handled by including the cross terms in the explicit part of the time splitting, 
leaving the Helmholtz problem for the velocity components uncoupled. 
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As seen from (29), the calculation of the VMS LES model terms requires sev- 
eral additional operations. The increase in total computational work will vary 
with the size and complexity of the simulation, but for the cases considered 
in this paper the increase is in the range 20-40%, with the smallest relative 
increase for the largest simulations. To put these numbers into perspective, we 
note that the total computational complexity of the spectral element method 
is OlK^N"^), so increasing the polynomial order {N — 1) from 6 to 7 gives a 
70% increase in computational time, about the same as increasing the number 
of elements in each dimension (K) from 5 to 6 would give. 



3.3.4 Smagorinsky model 

The eddy viscosity VT{x,t) is chosen in [1] as a Smagorinsky-type function: 

vt = {C'sMf\V'ul (30) 

or alternatively 

VT^{C'sMf\V'u\. (31) 

The former was labeled "small-small" in [2] , while the latter was labeled "large- 
small" . 

As the purpose of the model term is to approximate the effect of the unresolved 
scales on the small scales, it is argued in [1] that (30) is more consistent with 

the physical basis of the method, whereas (31) appears to be a computationally 
attractive alternative. The results in [2,3] show that good results are obtained 
with both methods. However, in terms of the spectral element implementation, 
the "large-small" form is not a computational simplification. A more attractive 
form is instead the "full-small" term 

i/T = (C^A')^|V^u|, (32) 

in which the scale extraction operators are avoided completely. 
The sum iV'^ttl can be written out as 



The constant C'g is set to 0.1, following [2,3], while A' is calculated for each 
element as the geometric average of the mean grid spacing in each direction. 
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4 Computational results 



4-1 Channel flow 

The plane turbulent channel flow is one of the simplest cases of an inhomoge- 
neous turbulence field, and this configuration has therefore been extensively 
used to assess the performance of turbulence models. The fully developed, sta- 
tistically stationary, plane channel flow is an equilibrium flow, because there 
is a global balance between the rate of production of turbulent kinetic energy 
and the rate of viscous dissipation. 

The fluid domain is bounded by two infinitely large parallel solid walls, and 

the flow is driven by a constant mean pressure gradient in the stream-wise di- 
rection along the walls. The boundary conditions are no-slip at the solid walls, 
and periodicity is imposed in the streamwise (x) and spanwise (z) directions, 
respectively. The wall-normal direction is thus y, and the channel half-height 
is denoted h. 

The instantaneous flow fleld is three-dimensional and time dependent, the 
ensemble (or time) averaged flow field is however unidirectional. If we let (•) 
denote the ensemble average, we therefore have U = (u) — [U{y), 0, 0]. 

The friction velocity, Ut, is defined by 



and this is used in the definition of the relevant Reynolds number for plane 
channel fiow: Re,- = Urh/u. 

Integrating the ensemble averaged Navier-Stokes equations in the wall-normal 
direction yields 



where the pressure gradient is a constant, related to the Reynolds number by 



Hence, the sum of the viscous {/idU/dy) and turbulent {—p{u'v')) stresses must 
vary linearly across the channel. The turbulent stress contribution dominates 
across the channel except very close to the wall where viscous stress dominates. 
This region is usually referred to as the viscous sub-layer and its thickness 
decreases with increasing Reynolds numbers. 




(35) 




(36) 
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Table 1 

Grid parameters for the present DNS and the reference simulations by Moser et 
al. [21] and by del Alamo et al. [22,23]. Grid spacing in wall units are calculated 
from the nominal Re,-. 
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Table 2 

Grid parameters for the VMS LES runs. Grid spacing in wall units are calculated 
from the nominal Re,-. 



We consider three different Reynolds numbers: Re,- = 180, 550, 950, and the 
VMS LES computations are compared with reference solutions obtained from 
direct numerical simulations. 
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Fig. 4. Details of the point and element distribution in the wall-normal direction 
for the grids listed in Tables 1 and 2. The longer bars show element boundaries for 
the spectral element grids. 

4-2 Direct numerical simulations at Re^ = 180 

As a first step towards our ultimate goal, to implement and evaluate the 
variational multiscale LES method in a high order spectral element fiow solver, 
we performed a Direct Numerical Simulation to verify the code. To this end, 
we considered fully developed channel flow at Re,- = 180, which corresponds 
to the well-known benchmark simulations reported by Kim et al. [24]. We 
performed the actual comparison of the results with the updated data set 
reported by Moser ct al. [21] who used a fully spectral Fourier/ Chebyshev 
method with 128 x 129 x 128 grid points. 

The simulation was carried out on a computational domain that approxi- 
mately corresponds to the one used by the reference solutions [21,24], see 
Table 1 for details. Across the channel we used 16 non-uniformly distributed 
elements with 8 nodal points in each element. In the streamwise and span- 
wise directions we used 16 x 16 uniformly distributed elements with 8x8 
nodal points per element. Thus, the total number of nodal points amounts 
to 112 X 113 X 112 in the streamwise, wall-normal, and spanwisc directions, 
respectively. The solution was advanced in time with a time-step correspond- 
ing to 0.18 viscous time-units and with 50% polynomial filtering [18] 
on each time-step. The simulation was initiated by a flow fleld obtained from 
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Fig. 5. Rct- = 180: Variation of the mean velocity across half the channel in viscous 
units, compared with the reference solution of Moser et al. [21]. 

an existing plane channel flow solution obtained by a finite-volume code. The 
flow then evolved approximately 54 flow-through times before a fully devel- 
oped state was achieved. The results presented here was obtained by collecting 
statistics over approx. 20 flow-through times. The flow statistics are averaged 
over the homogeneous - strcamwise and spanwise - directions. Homogeneity 
in a speciflc direction implies that any correlation of a fluctuating quantities 
remains invariant under translation in that direction. 



4.2.1 Results 

The actual computed Reynolds number is Re* = 178.83, i.e. within 0.7% of 
the prescribed value and well within what can be expected. Moser et al. [21] 
reported Re* = 178.13. The results presented in Figs. 5-8 compare in all 
aspects very well with the benchmark data, thus estabhshing solid confldence 
in the numerical method. The slight deviations reported herein is well within 
what should be expected, and even closer correspondence could have been 
obtained by simply collecting statistics for a longer period of time. This was, 
however, not considered to be necessary. 

As background for the VMS LES results presented below, we also include 
results from a simulation on the grid "Coarse-36" (see Table 2 for grid prop- 
erties). This simulation contains no turbulence modelling, but 2% polynomial 
flltering [18] is employed. Except for the pressure correlations in Fig. 8, the re- 
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Fig. 6. Re,- = 180: Variation of mean viscous shear and the turbulent shear stress 
across half the channel, compared with the reference solution of Moser et al. [21]. 
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Fig. 7. Re,- = 180: Variation of streamwise {u'), spanwise {v'), and wall-normal {w') 
root-mean-square velocity fluctuations across half the channel, compared with the 
reference solution of Moser et al. [21]. 
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Fig. 8. Re,- = 180: Variation of the root-mean-square pressure fluctuations across 
half the channel, compared with the reference solution of Moser et al. [21]. 

suits are so good that modelling is not expected to improve them. This shows 
that the spectral element method gives high accuracy even for relatively coarse 
grids, but it also indicates that plane channel flow is not the most challenging 
test case. The availability of quality reference data makes it attractive as a 
starting case, we must however keep in mind that the grids for the model tests 
have to be sufficiently coarse and not turn into a "quasi-DNS" e.g. near the 
walls. 



4.3 VMS LES results 



Lots of combinations of the scale partitioning parameter and the Smagorinsky 
forms were tested for Re^ = 180, and the best choice was used for additional 
simulations at Re^ = 550 and Re^ = 950. 

The spectral element grid for Re^ = 180 was chosen as the "Coarse-24" grid 
described in Table 2. The element interfaces in the wall-normal direction were 
given by a coarse Gauss-Lobatto-Chebyshev grid, as recommended in [25]. 
The scale partitioning cut-off mode was kept constant for all elements, even 
though the element size varied in the wall-normal direction. 

In order to get a real test of the modelling, the grid had to be much coarser than 
what would give reasonably good results without a model. Spectral element 
grids for the higher Reynolds numbers were constructed such that the first 
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Fig. 9. Rct- = 180: Variation of the mean velocity across half the channel in viscous 
units, compared with the reference solution of Moser et al. [21]. 

element interface in the wall-normal direction is placed at approximately the 
same value of for all the cases, see the illustrations in Fig. 4. To reduce the 
number of parameters, the polynomial degree was fixed for all the VMS LES 
runs; only the number of elements was changed. 



4.. 3.1 Simulations at Re^ — 180 

The grid parameters for this case are given in the column "Coarse-24" in 
Table 2. 

Without a model, both over-integration and polynomial filtering (2%) was 
necessary to keep the simulation stable at this resolution. With the VMS model 
term, either method was sufficient. It was found that polynomial filtering did 
reduce rather than improve the quahty of the results. To obtain the presented 
VMS results we therefore employed only over-integration in the simulations. 

Beside using the different forms of the Smagorinsky term (30)-(32), the scale 
partitioning was varied in the simulations. With a local grid of 7 grid points in 
each direction on each element, we have used N — A and N — h ior the large- 
scale extraction described in Section 3.3.1. These values correspond to 57% and 
71% of the one- dimensional spectrum, respectively. In three dimensions, the 
resulting large-scale spaces consist of 19% and 35% of the modes, respectively. 

Varying the scale partitioning had a strong influence on the results, and N — 5 
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Fig. 10. Re,- = 180: Variation of streamwise («'), spanwise {w'), and wall-normal 
{v') root-mean-square velocity fluctuations across half the channel, compared with 
the reference solution of Moser et al. [21]. 
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Fig. 11. Rcr = 180: Variation of the mean velocity across half the channel in viscous 
units, compared with the reference solution of Moser et al. [21]. 

was found to be the best choice, as seen from Figs. 9 and 10. The rest of the 
results shown here are obtained with N — 5. 
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Fig. 12. Re^ = 180: Variation of mean viscous shear and the turbulent shear stress 
across half the channel, compared with the reference solution of Moser et al. [21]. 
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Fig. 13. Re^ = 180: Variation of streamwise {u'), spanwise {w'), and wall-normal 
(v') root-mean-square velocity fluctuations across half the channel, compared with 
the reference solution of Moser et al. [21]. 
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Fig. 14. Re-r = 180: Variation of the root-mean-square pressure fluctuations across 
half the channel, compared with the reference solution of Moser et al. [21]. 
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The different forms of the Smagorinsky term gave very similar results for Re^ — 
180. The results are presented in Figs. 11-14. The results from "large-small" 
form (31) were almost indistinguishable from the "full-small" (32) results, and 
are not included in the figures. 

As shown in Section 4.2.1, simulations on the "Coarse-36"-grid gave good re- 
sults without modelling for this case. Results from simulations without mod- 
elling on an intermediate grid with 30'^ grid points were comparable to the 
VMS results from the 24^-grid shown here, but at a 40% higher computa- 
tional cost. 



Simulations at Re,- = 550 

The grid parameters for this case are given in the column "Coarse-42" in 
Table 2. 

The scale partitioning parameter of = 5, which was found to be the best 
choice for Re,- = 180, was also used for this case. Again, the "full-small" and 
"large-small" Smagorinsky forms produced very similar results, so the latter 
are not shown. The results are presented in Figs. 15-17. 
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Fig. 15. Re,- = 550: Variation of the mean velocity across half the channel in viscous 
units, compared with the reference solution of del Alamo and Jimenez [22]. 
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Fig. 16. Rbt- = 550: Variation of mean viscous shear and the turbulent shear stress 
across half the channel, compared with the reference solution of del Alamo and 
Jimenez [22]. 
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Fig. 17. Re,- = 550: Variation of streamwise {u'), spanwise {w'), and wall-normal 
{v') root-mean-square velocity fluctuations across half the channel, compared with 
the reference solution of del Alamo and Jimenez [22]. 
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Fig. 18. Rc,- = 950: Variation of the mean velocity across half the channel in viscous 
units, compared with the reference solution of del Alamo et al. [23]. 
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Fig. 19. Re^ = 950: Variation of mean viscous shear and the turbulent shear stress 
across half the channel, compared with the reference solution of del Alamo et al. [23]. 
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Fig. 20. = 950: Variation of streamwise {u'), spanwise {w'), and wall-normal 
(v') root-mean-square velocity fluctuations across half the channel, compared with 
the reference solution of del Alamo et al. [23]. 
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4-3.3 Simulations at Re,- = 950 

The grid parameters for this case are given in the column "Coarse-60" in 
Table 2. 

In this case we have only run the "full-small" Smagorinsky form, and the scale 
partitioning parameter is still N — 5. The reference simulation is described 
in [23], and the reference data arc downloaded from the site given in [22]. Our 
results are presented in Figs. 18-20. 



4-4 Comments on the results 



The VMS LES results show clear improvement from the results without a 
model, in particular for the higher Reynolds numbers. The plane channel flow 
at Re^ = 180 does not seem to provide sufficient challenges for the testing 
of turbulence models, as it is too easy to resolve the main features without 
any modelling at all. The VMS LES results are not compared with alterna- 
tive turbulence models, as the intentions of this study was mainly to lay the 
foundations for the incorporation of VMS LES in a spectral element method. 
Therefore only the simplest Smagorinsky eddy viscosity was used in the model 
terms in the small-scale equations. 



5 Conclusions 



The Variational Multiscale Large Eddy Simulation method has been imple- 
mented within the framework of a spectral element method. The presented 
scale partitioning method was shown to produce a gradual introduction of the 
small-scale model terms. This is intuitively favourable to a sharp cut-off at a 
given point in the spectral space. The computational overhead for the method 
was 20-40% for the applications considered here. This must be considered 
to be reasonably low, as even small increases in the spatial resolution of the 
spectral element method are more computationally demanding. Good results 
have been obtained for plane channel flows up to Re,- = 950, even for grid 
densities as low as 0.06% of the reference simulation grid density, and using 
the simplest possible small-scale dissipation model. 
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